// P21ccDiffractor.cpp: implementation of the P21ccDiffractor class.
//
//////////////////////////////////////////////////////////////////////


#include "P21ccDiffractor.h"

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

///////////////////////////////////
// b unique
///////////////////////////////////
P21ccDiffractor::P21ccDiffractor(char *DiffractorFile, 
						   Position *W1, Position *W2, Position *W3, char **atomLabel, int nAtom)
: Diffractor(DiffractorFile, W1, W2, W3, atomLabel, nAtom)
{
newA(Case, NIntensity);
newA(a1,   NIntensity);
newA(a2,   NIntensity);
newA(a3,   NIntensity);
newA(AP,   NAtom);

for(int n=1;n<=NIntensity;n++)
	Case[n] = ( (K[n]+L[n])%2 == 0 );
}


P21ccDiffractor::~P21ccDiffractor()
{
deleteA(Case);
deleteA(a1);
deleteA(a2);
deleteA(a3);
deleteA(AP);
}


//////////////////////////////////////////////////////////////////////
// functions
//////////////////////////////////////////////////////////////////////

void P21ccDiffractor::TransformDerivatives(double   *IDerivative, 
										   Position *XDerivative)
{
int    n, m, h, k, l;
double x, p, q, F1, F2, F3, *FF;

for(n=1;n<=NIntensity;n++)
	{
	x = IDerivative[n]*ReAmplitude[n];

	a1[n] = DH[n]*x;
	a2[n] = DK[n]*x;
	a3[n] = DL[n]*x;
	}

for(m=1;m<=NAtom;m++)
	{
	SetUpTrig(AP[m]);

	FF = FiniteFactor[m];

	F1 = 0.0;
	F2 = 0.0;
	F3 = 0.0;

	for(n=1;n<=NIntensity;n++)
		{
		h = H[n];
		k = K[n];
		l = L[n];

		p = FF[n]*CosY[k]*(SinX[h]*CosZ[l] + CosX[h]*SinZ[l]);
		q = FF[n]*SinY[k]*(CosX[h]*CosZ[l] - SinX[h]*SinZ[l]);

		if(Case[n])
			{
			F1 += a1[n]*p;
			F2 += a2[n]*q;
			F3 += a3[n]*p;
			}

		else
			{
			F1 += a1[n]*q;
			F2 += a2[n]*p;
			F3 += a3[n]*q;
			}
		}

	FDeriv[m][1] = -F1;
	FDeriv[m][2] = -F2;
	FDeriv[m][3] = -F3;
	}
			

MakeXDerivative(XDerivative);	
}



void P21ccDiffractor::FourierTransform(Position *AtomPosition)
{
double *FF;

int    n, m, h, k, l;

for(m=1;m<=NAtom;m++)
	AP[m] = AtomPosition[m];


for(n=1;n<=NIntensity;n++)
	{
	ReAmplitude[n] = 0.0;	
	ImAmplitude[n] = 0.0;
	}


for(m=1;m<=NAtom;m++)
	{
	SetUpTrig(AtomPosition[m]);

	FF = FiniteFactor[m];

	for(n=1;n<=NIntensity;n++)
		{
		h = H[n];
		k = K[n];
		l = L[n];

		if(Case[n])
			ReAmplitude[n] += FF[n]*(CosX[h]*CosZ[l]-SinX[h]*SinZ[l])*CosY[k];
		else
			ReAmplitude[n] -= FF[n]*(SinX[h]*CosZ[l]+CosX[h]*SinZ[l])*SinY[k];
		}
	}
}